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ABSTRACT 

We study the evolution of the halo-halo correlation function and bias in four cosmological models 
(ACDM, OCDM, tCDM, and SCDM) using very high-resolution iV-body simulations with dynamical 
range of ~ 10,000 — 32,000 (force resolution of w 2 — 4/i _1 kpc and particle mass of sa 10 9 h~ 1 M Q ). 
The high force and mass resolution allows dark matter (DM) halos to survive in the tidal fields of high- 
density regions and thus prevents the ambiguities related with the "overmerging problem." This allows 
us to estimate for the first time the evolution of the correlation function and bias at small (down to 
~ 100/i _1 kpc) scales. 

We find that at all epochs the 2-point correlation function of galaxy-size halos £hh is well approximated 
by a power-law with slope w 1.6 — 1.8. The difference between the shape of £hh and the shape of the 
correlation function of matter results in the scale- dependent bias at scales <J 7h~ 1 Mpc, which we find to 
be a generic prediction of the hierarchical models, independent of the epoch and of the model details. 
The bias evolves rapidly from a high value of ^2 — 5 at z~3 — 7 to the anti-bias of b ~ 0.5 — 1 
at small <; 5/i _1 Mpc scales at z = 0. Another generic prediction is that the comoving amplitude of 
the correlation function for halos above a certain mass evolves non-monotonically: it decreases from an 
initially high value at z ~ 3 — 7, and very slowly increases at z <^ 1. We find that our results agree 
well with existing clustering data at different redshifts, indicating the general success of the hierarchical 
models of structure formation in which galaxies form inside the host DM halos. Particularly, we find an 
excellent agreement in both slope and the amplitude between £hh(z = 0) in our ACDM 6 o simulation and 
the galaxy correlation function measured using the APM galaxy survey. At high redshifts, the observed 
clustering of the Lyman-break galaxies is also well reproduced by the models. We find good agreement at 
z > 2 between our results and predictions of the analytical models of bias evolution. This indicates that 
we have a solid understanding of the nature of the bias and of the processes that drive its evolution at 
these epochs. We argue, however, that at lower redshifts the evolution of the bias is driven by dynamical 
processes inside the nonlinear high-density regions such as galaxy clusters and groups. These processes 
do not depend on cosmology and tend to erase the differences in clustering properties of halos that exist 
between cosmological models at high z. 

Subject headings: cosmology: theory - large-scale structure of universe - methods: numerical 



1. INTRODUCTION 

It is widely believed that the distribution of galaxies 
is different from the overall distribution of dark matter 
(DM). This difference, the bias, is crucial for comparisons 
between observations and predictions of cosmological mod- 
els. Observations provide information about the distribu- 
tion of objects such as galaxies and galaxy clusters. The 
models, however, most readily predict the evolution of 



the dark matter diotribution, which cannot bo oboorvod 
directly. The models, therefore, should be able to pre- 
dict the distribution of objects or, conversely, predict how 
this distribution is different from that of the dark mat- 
ter (i.e., the b ias). The notion of bias was introduced by 



Kaiser (1984) to explain the large difference in cluster- 
ing strength b etw een galaxies and Abell clusters. Later 
Kaiser (1986) and Bardeen (1986) applied this argument 
to galaxies themselves. Davis et al. (1985) showed that 



the CDM model disagreed with observations, if galaxies 
were distributed like dark matter. However, if a biased 
galaxy formation scenario was assumed, the "galaxy" cor- 
relations substantially exceeded the correlations of mass at 
all scales and agreed with observations for certain values 
of bias. This work was followed by other studies which 
tried to account for the phenomenon of gal a xy bias (e.g., 
iRocs 19851 ; |Schacffer fc Silk 1985fc |Silk 1985| ; |Dckel fc Silk 



1986) 



The bias can be defined and understood differently. In 
this paper we will use the conventional statistical defini- 
tion of the bias as the ratio of the correlation functions of 
objects and dark matter: 



b 2 (M, r,z) 



= tgg( M , r, z) 
idm{r, z) 



(1) 



Here, b 2 is the square of the bias function, £ sg (M, r, z) and 
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£dm {r, z) are the 2-point spatial correlation functions of ob- 
jects and dark matter, respectively. Dependencies in the 
above equation indicate that in general the bias may de- 
pend on the epoch z, scale r, and properties of the objects 
such as their mass M. It is also expected that the func- 
tional form b 2 (M, r, z) depends on the cosmological model. 
The bias in the above definition is most closely related to 
the observations because the observed ^ gg (M, r, z) can be 
compared to the £d m (r, z) predicted within a framework of 
a given cosmological model. This gives an estimate of the 
amount of bias needed for a cosmological model to agree 
with observations. While it is clear that this approach 
cannot be used to test a model, it provides some insights 
into the nature of bias and its evolution. For example, 
Steidel et al. (1998) used the observed clustering strength 
of high-redshift (z k, 3) Lyman-break galaxies to derive 
the implied value of bias in different cosmological mod- 
els. This value was found to be quite large in all models 
(b ~ 3 — 6). On the other hand, a number of theoreti- 



cal stidies (e.g., Klypin, Primack & Holtzman 1996; Cole 
et al. 1997; Jenkins et al. 1998; and references therein) 
estimated t^dm in different cosmological models at z = 
and made comparisons with accurate measurements of £ gg 
from local galaxy surveys. These studies indicate that sig- 
nificant anti-bias (b < 1) is required for the open and flat 
low-f2o models at scales 3 — 8/i _1 Mpc, while Oo = 1.0 
cluster- normalized models indicate positive bias (6 > 1). 
Comparison of the low- and high-z results implies that, 
regardless of cosmological model, the bias has decreased 
significantly from the early epochs to the present. 

The distribution of the dark matter and £d m (r, z) cannot 
be observed directly. Therefore, a test of a cosmological 
model is possible only if the model can predict the £ sg of 
observed objects. Unfortunately, this is not an easy and 
straightforward task. First of all, we should fully under- 
stand what are the observed objects and where/how they 
form. The standard lore is that observed galaxies form 
dissipatively inside dark matter halos. The properties of 
galaxies will then depend on the mass of the parent halo, 
its spin, details of dissipative processes and mass accre- 
tion history, and other factors (e.g., Mo, Mao & White 
1998). Therefore, in order to predict the type of galaxy 
and its properties, the relevant processes must be included 
in the model. However, it seems likely that in every suf- 
ficiently massive (M ;> 1O 11 /i _1 M0) gravitationally bound 
halo baryons will cool, form stars, and produce an object 
ressembling a galaxy (e.g., Kauffman, Nusser & St einmetz 
1997; Roukem a et al. 1997; Yepes et al. 1997; [Salucc^ 
& Peisic 1997). The galaxy population as a whole can be 



evolution of bias in such models can be calculated quickly 
which makes extensive parameter-space studies possible. 
The main disadvantages are large uncertainties (especially 
at small, r <J 5h~ 1 Mpc, scales) introduced by bias pre- 
scriptions and limited applicability of the non-linear clus- 
tering approximations. 

Direct numerical simulations should, ideally, predict 
halo-halo clustering without any additional assumptions 
and uncertainties. However, until very recently the predic- 
tions of numerical simulations were also quite uncertain. 
The main reason for the uncertainty was that dissipation- 
less A-body simulations had been consistently failing to 
produce galaxy-size dark matter halos in dense environ- 
ments typical for galaxy groups and clusters. Recently, 
it was shown that this effect, known as "the overmerg- 
ing problem" (e.g., Frenk et al. 198~8| ; Summers, Davis, 
& Evrard 1995), is due mainly to the insufficient force 
and mass r esolution of such sim ulations (Moore, Katz fc 
Lake 1996; |Klypin et al. 1998j hereafter KGKK; phigna 
1998| ). The lack of sufficient resolution leads to ar- 
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then viewed as a population of galaxy-size dark matter ha- 
los. The clustering of the latter can be studied without in- 
clusion of complicated physics; it has been modelled using 
both direct numerical simulations and analytical methods. 
Typical ingredients of the analytical models (e.g. 



Mata rrese e t al 
et al 



1997; Moscardini et al. 1998; 



Manrj 



1998) are the extended Press-Schechter formalism 



( Bower 1991 ; Bond ct al. 1991 ) used to follow mass evolu- 
tion of halos, analytical approximations to the non-linear 
clustering evolution of the dark matter (e.g., Hamilton et 
al. 1991; Peacock & Dodds 1994,1996; Jain, Mo & White 
1995; Smith et al. 1998), and a model for bias evolution 
(e.g., Mo & White 1996; Matarrese et al. 1997). The 



tificial disruption of halos in clusters. This, in turn, leads 
to a strong artificial anti-bias (especially at small scales 
<J 3/i _1 Mpc, but larger scales are also affected). There 
are several ways to deal with this problem. One possible 
way is to break up massive structureless halos into subha- 
los using s ome kind of observationally motivated prescrip- 
tion (e.g., Nolthcnius, Klypin, & Primack 1997; Klypin, 
Nolthcnius, & Primack 1997). These subhalos can then 
be included into halo catalogs used to compute the cor- 
relation function and other halo statistics. Another com- 
mon approach is to overcome overmerging by weighting the 
massive halos according to their ma ss and comp ute thus a 



weighted correlation function (e.g., Bagla 1998). Both of 
these approaches are useful. Nevertheless, it is not clear 
whether the unavoidable heuristic assumptions take the 
processes driving the small-scale bias evolution correctly 
into account. The weighting technique, for example, ig- 
nores the real physical effects in groups and clusters such 
as tidal stripping and dynamical friction. We will argue 
below that these effects are likely to be driving the small- 
scale bias evolution at low (z <; 1) redshifts. Hydrody- 
namic simulations that include gas cooling are affected by 
overmergin g to a significantly lesser degree (e.g ., Summers 
etal. 1995; [Katz, Hernquist, fe Weinberg 199S| ). The cool- 
ing creates compact dense objects inside halos which can 
survive in clusters. These simulations, therefore, can be 
used to study the halo clustering directly. Unfortunately, 
the computational cost required to simulate a large vol- 
ume with sufficiently high mass resolution is prohibitively 
high. Moreover, such simulations usually oversimplify the 
gas dynamic by including only cooling mechanism. This 
leads to "overcooling" : without a heating process to regu- 
late it, the cooling produces very compact and dense bary- 
onic blobs in the halo centers. These blobs do survive suc- 
cessfully in clusters, but they also suffer much less from 
the tidal stripping of material as compared to a realistic 
galaxy with a more extended distribution of baryons. The 
mass of the objects may thus be higher than it should be 
which may lead to excessive dynamical friction and thus 
incorrect dynamics of halos. 

The very high resolution A-body simulations are thus a 
viable alternative, if the required resolution can be reached 
at an affordable computational cost. Analytical arguments 
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and numerical experiments (Moore et al. 1996; KGKK) 
indicate that the required force and mass resolution are 
s=a 1 — 2h~ 1 kpc and <; 10 9 /i _1 M Q , respectively. This force 
resolution is sufficiently high for most halos to survive in 
high-density regions. The mass resolution is determined 
by the requirement that a galaxy-size halo should have 
;> 100 particles to avoid relaxation effects and to assure 
a robust identification by a halo finding algorithm. The 
numerical simulations that reach such resolution (Ghigna 
et al. 1998; KGKK) show that most halos do survive even 
in the richest clusters. The dynamic range required to 
reach the needed resolution in a statistically large volume, 
- 50 - lOO/i^Mpc, is quite high: 2 - 5 x 10 4 . Never- 
theless, the advances in computer hardware and in the 
numerical algorithms make it now possible to carry out 
such simulations at an affordable computational cost. In 
this paper we use such high dynamic range simulations 
to calculate the 2-point correlation function of halos and 
corresponding bias using all halos in the simulated vol- 
ume: i.e., both isolated dark matter halos and satellites 
of massive halos, and halos inside group- and cluster-size 
systems. The simulations of different cosmological models 
and box sizes were made using the Adaptive Refinement 
Tree (ART; Kravtsov, Klypin & Khokhlov 1997) and the 
AP 3 M (Couchman 1991) A-body codes. The absence of 
the overmerging 1 in these simulations means that the ad- 
ditional steps such as breaking-up of clusters, or mass- 
weighting (see above) are not necessary. Therefore, the 
correlation function of halos is measured directly down to 
unprecedentedly small scales (w 150/i _1 kpc) without the 
usual uncertainties associated with these steps. For the 
first time this opens the possibility to study the effects 
of dynamical processes such as tidal destruction and dy- 
namical friction on the evolution of small-scale bias. The 
results on the evolution of bias and on the halo correla- 
tion function can be used as a basis for comparisons and 
interpretations of the existing and upcoming observations, 
as well as a check and/or input for the analytic models of 
clustering evolution. 

The paper is organized as follows. In § 2 we briefly re- 
view the definitions of bias and current analytical models 
of its evolution. The cosmological models studied in this 
paper are described in § 3. The details of the numerical 
simulations and discussion of the construction of halo cat- 
alogs and halo survival in the high-density regions is given 
in § 4. In § 5 we present our results on the evolution of 
halo clustering and bias in different cosmological models. 
We also present a comparison of our z = results with 
the galaxy correlation function measured using the APM 
galaxy survey. A discussion of the main results is presented 
in § 6. We summarize our main results and conclusions in 
§ 7. 

2. THE NOTION OF BIAS 

The notion of bias is more complicated than eq. (1) 
might suggest. Formally, the bias is defined as a function 
relating fluctuations in the dark matter density, 8d m , to 

The extent to which the very central regions (r < 200/i _1 kpc) of clusters are affected by the overmerging, even with resolution this high, 
is still a matter of debate. Ghigna et al. (1998), for example, argue that these regions are still completely overmcrged. They find, however, 
that halos survive at smaller radii (as 50/i~ 1 kpc) in their 5 kpc resolution RUN1, as compared to the 10 kpc RUN2. The resolution of the 
simulations presented in this paper is comparable to the RUN1. We do find halos within central 100/i _1 kpc (see § 4.3 and Figs. 1-3). 
Although some halos do survive at r < 100ft -1 kpc from cluster center, our tests show that these regions may be affected by the overmerging. 

The larger scales > lOO/i kpc, are not affected 



the fluctuations in the number density of objects, S n . In 
general, this relationship can be complicated and may de- 
pend on a large number of factors: scale, mass and/or type 
of objects, time, etc. The large number of the dependen- 
cies or discreteness of the object distribution can result 
in the stochasticity of the bias (Dekel & Lahav 1998): a 
scatter in the relationship between S n and 5 dm- Finally, 
the bias can be non-local, if probability to form an object 
at a given point is not fully determined by local factors. 
For example, in addition to the local density, temperature, 
etc., the probability to form a galaxy may depend on en- 
vironment. Analytical results of Catelan et al. (1998b) 
indicate that the bias is expected to be non-local even in 
the linear regime. 

The lack of general understanding of these dependen- 
cies of the bias usually results in the use of the simplest 
assumptions. The most common approach is to assume 
that bias is local, depends only on the local matter den- 
sity, and is linear. S n = bSdm- An obvious consequence 
of the latter assumption is that bias is scale- independent. 
In this case, the correlation function and power spectrum 
of halos and DM are simply related by a constant scaling 
factor. Although the linearity of bias would greatly sim- 
plify the theoretical interpretation of the clustering data, 
it is likely that bias depends on a variety of processes that 
may lead to nonlinearity. These processes include merg- 
ing, tidal disruption (e.g., Dubinski 1998), suppression of 
galaxy formation in small halos due to supernovae feed- 
back (Dekel & Silk 1986; Yepes et al. 1997), etc. The 
observations, in fact, indicate that at small scales galac- 
tic bias is nonlinear. The correlation functions of different 
types of galaxies differ in amplitude and shape suggesting 
that at least some of the galaxies are nonlincarly biased. 

During the last years, a significant progress has been 
made in the analytical modelling of the bias and its evo- 
lution (e.g., Coles 1993; Fry 1996; Mo & White 1996, 
hereafter MW; Matarrese et al. 1997; Mann et al. 1997; 
Catelan et al. 1998a, b). Although the current analytical 
models may have some drawbacks and limitations, they 
provide an insight and interpretation for numerical simu- 
lations. At this point it is also important to check how 
well the prediction of the analytical models agree with the 
results of simulations (see, e.g., MW; Mo, Jing & White 
1996; Jing 1998). We will therefore review briefly some re- 
sults of the analytical models concerning evolution of the 
bias. 

In a seminal paper, Kaiser (1984) showed that if the 
observed systems in the universe (such as galaxies and 
galaxy clusters) form in the peaks of the density field, their 
distribution is biased. This is an example of statistical 
bias: the bias determined simply by the fact that objects 
do not sample the distribution of matter but that of the 
peaks, the latter having a statistically different distribu- 
tion from the former. Thus, the bias is introduced from 
the start by the way the objects form. However, the bias 
changes subsequently as the distribution evolves driven by 
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the gravitation. In hierarchical models this evolution re- 
sults in growth of objects via multiple mergers. However, 
if the merger rate is low, the evolution of bias can be mod- 
eled by the object- conserving model. In this model the 
objects form with an initial statistical bias and after that 
are dragged without merging by a gravitational pull from 
the surrounding density fluctuations (Dekel & Rees 1987; 
Nusser & Davis 1994; Fry 1996). 

MW used the extended Press-Schechter formalism to de- 
rive an expression for the bias in Lagrangian coordinates 2 
(comoving radius R of the region from which halos form 
or halo mass M). At linear scales, this expression 3 for 
the bias at redshift z for DM halos of mass M formed at 
redshift Zf is (Matarrese et al. 1997): 



b(M,z\z f ) = 1 + 



1 



(2) 



Here 5f — S c D + (z)/D + (zf), D + (z) is the growing mode 
of linear perturbations normalized to unity at z = 0, 
v = Sf/a(M,z), 5 C is the critical overdensity for spher- 
ical collapse at z = 0, and <r(M, z) is the rms linear mass 
fluctuation on the scale M of halos linearly extrapolated 
to redshift z. This expression is similar to that obtained 
in earlier studies from a peak-background split argument 
(e.g., Cole & Kaiser 1989). The standard interpretation 
of the Press-Schechter description of the hierarchical evo- 
lution is that at any epoch z all halos merge immediately 
to form more massive halos. Thus, if observed objects 
are identified with host halos at any epoch, then z — Zf 
in eq. (2). Matarrese et al. (1997) call this the merging 
model. The object-conserving and merging models are two 
extremes pictures of clustering evolution, although they 
may be applicable to the evolution of galaxy clustering 
at certain epochs. At some epochs the halos may neither 
survive nor merge instantly. Moreover, other processes 
such as halo dynamics in galaxy clusters and groups may 
become important when clustering reaches highly nonlin- 
ear stages. The wealth of potentially important processes 
may make a one-to-one identification between galaxies 4 
and DM halos very difficult (see, for example, discussion 
in Moscardini et al. 1998). 

Eq. (2) gives an estimate of the bias for objects of a 
single mass. The bias of a sample of objects ("effective" 
bias 5 ) with a range of masses M > M min should be calcu- 
lated as a weighted average over the mass distribution of 
objects n(M, z): 



b{M,z\zf)n(M,z)dhxM, (3) 



M>M„ 



where n(M, z) is given by the Press-Schechter (1974) dis- 
tribution and n(z) is the mean number density of objects 
with masses M > M m i n at redshift z. Moscardini et al. 
(1998) give a useful fitting formula for b e ff(z): 



b eff{z) 



l/S c +[b eff (Q)-l + l/5 c ]/D + (zr, (4) 



and provide best-fit parameters 6 e //(0) and /? for a vari- 
ety of different cosmological models and values of M min . 
The generic feature of the evolution of the effective bias 
described by eq. (4) is its rapid decrease with decreasing 
redshift and increase with increasing M m i n . For galaxy- 
size halos, b eff (0) ~ 0.5 - 1.0 and /3 ~ 2.0 - 1.7 for 
M min ~ 10 9 - lO 12 ft. _1 M . This gives b eff ~ 2 - 4 at 
z = 3 and b e ft ~ 0.5 — 1.0 at the present epoch. We will 
use eq. (4) in § 6 to interpret and compare our results 
with predictions of the analytical models discussed above. 



3. COSMOLOGICAL MODELS 

We have chosen to study the evolution of bias in four 
representative variants of the CDM models (see details 
in Table 1): 1) the standard COBE-normalized CDM 
model (SCDM); 2) a variant of the fi = 1.0 CDM model 
(rCDM) with a different shape of the power spectrum (the 
shape parameter T = floh = 0.2); 3) a flat low-density 
model with fi Q = 1-Q\ = 0.3 (ACDM); 4) an open model 
with f2o = 0.3 (OCDM). The observations of the galaxy 
clustering indicate that the power spectrum of the galaxy 
distribution has a shape different from that of the stan- 
dard CDM model: T = n Q h « 0.2 instead of T = 0.5 of 
the SCDM (e.g., Maddox, Efstathiou, & Sutherland 1996; 
Peacock & Dodds 1994). This has motivated Jenkins et 
al. (1998) to study the rCDM model, in which r indicates 
the fact that lower values of T may be obtained with a 
late decay of the massive T-neutrino. The physics (addi- 
tional to the physics of the SCDM model) responsible for 
the change of the shape is, in fact, irrelevant for study of 
structure formation. Although the model in that respect 
is somewhat heuristic, it is interesting as an example of a 
model with flo — 1 and an approximately correct shape 
of the power spectrum. The observations of the galaxy 
cluster evolution (Eke et al. 19 98) and of the baryon frac- 
tion in clusters ( Evrard 1997| ) strongly indicate value of 
matter density f2o ~ 0.3, while various observational mea- 
surements of the Hubble constant (e.g., Kim et al. 1997; 
Falco et al. 1997; Salaris & Cassisi 1998) tend to con- 
verge on the values of ft « 0.6 — 0.7. Therefore, we have 
considered two models (open and flat) with = 0.3 and 
ft = 0.65 — 0.7. The age of the Universe in all of our models 
is given in Table 1 and is in good agreement with the ages 
of the oldest globular clusters (Chaboyer 1998). 

We have used different approximations for the power 
spectrum of density fluctuations for different considered 
models. For the OCDM and ACDMi models (see Table 1) 
we used the Bardeen et al. (1986; BBKS) fit for the power 
spectru m P(k) — AkT 2 (k) with corrections of Sugiyama 



(1995 



T{k) 



ln(l + 2.34q) 

2.34g 
(I6.L3) 2 



(5.46 g ) E 



3.89g + 
f (6.7L3) 4 ]- 1 / 4 



(5) 



2 More sophisticated analytical treatment of bias in Lagrangian coordinates can be found in Catelan et al. (1998b) and Mann et al. (1997). 

3 MW note that equation (2) should be valid even when <5<j m > 1 or £,dm( T ) ~ 1, as long as the scale r is larger than the Lagrangian radius 
R. They find that at these scales the correlation function of halos and the bias in their numerical simulations are well described by equation 
(2) (see, however, Catelan et al. 1998a, and Jing 1998). 

4 The situation with evolution of clusters of galaxies is, of course, much simpler. 
5 We follow here the notation and terminology of Matarrese et al. (1997). 
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Table 1 
cosmological models 



Model 


^0 


^A,0 


h 


to 


08 


Approximation 








(Gyr) 






SCDM 


1.0 


0.0 


0.50 


13.1 


1.1 


Efstathiou et al. 


OCDM 


0.3 


0.0 


0.65 


12.2 


0.9 


BBKS+Sugiyama 


ACDMi 


0.3 


0.7 


0.70 


13.4 


1.17 


BBKS+Sugiyama 


ACDM 2 


0.3 


0.7 


0.70 


13.4 


1.0 


Klypin & Holtzman 


rCDM 


1.0 


0.0 


0.5 


13.1 


1.0 


Efstathiou et al. 



and 



Qoh" exp[— fib — V2h(Qi,/tto 



(6) 



where £l b = 0.0125/i~ 2 flWalker et al. 1991p ). The power 



spectra for both the SCDM and the r CDM models were 



appro ximated by the fitting formula of Efstathiou, Bond 



& Wl litc (1992) 



P(k) 



Ak 



1 + [ak + (bkf 2 + (ck) 2 ] 



211.13 N l 2 / 113 



(7) 



where a = 6.4/r, b = 3.0/r, c = 1.7/T, and A is the nor- 
malization constant. The shape parameter T is 0.5 and 
0.2 for the SCDM and tCDM models, respectively. These 
two analytic fits provide fairly good approximations to the 
power spectra of these models in the limit f2b/i7o -C 1- 

For the ACDM2 models we have used an approximation 
to the power spectrum different from that of the ACDMi 
model. The approximation 



P(k) 



Ak 



(l + axk 1 ' 2 + a 2 k + a 3 k 3 / 2 + a 4 k 2 ))' 



(8) 



where ai = -1.5598, a 2 = 47.986, a 3 = 117.77, a 4 = 
321.92, and CI5 = 0.9303, is given by Klypin & Holtzman 
(1997) and was obtained by a direct fit to the power spec- 
trum estimated using a Boltzmann code. The accuracy 
of this approximation is < 2% (see Klypin & Holtzman 
1997 for details). A small shift in the normalization makes 
the approximations for the ACDMi and ACDM2 models 
very similar in the range of wavenumbers probed in our 
simulations (k w (0.1 — 10)/iMpc _1 ). This shift can be ex- 
pressed using the value of as (amplitude of fluctuations on 
8/i _1 Mpc scale): as = 1-17 for the BBKS approximation 
versus as — 1-0 for the Holtzman approximation. With 
these values of as the differences in power are negligible at 
large scales (k « (0.1 — 0.3)/iMpc~ ) and are within 10% 
at the smaller scales. 

Our SCDM model was normalized to the two-year 
COBE-DMR data, as = 1.1 (a somewhat higher value, 
as ~ 1.15 — 1.2, is obtained from the four- year COBE- 
DMR data, see e.g., Bunn & White 1997). This normal- 
ization is inconsistent with normalization, as ~ 0.5 , de- 
duced from the observed cluster abundances (e.g., Eke, 

6 For the ACDM2 we used a slightly higher value of Q b = 0.015?t -2 . 



Cole, fc Frenk 1996] ) , which reflects the well-known failure 



of this model to account for both COBE and cluster data. 
Our normalization of the OCDM is higher than that im- 
plied by the four-year COBE-DMR data (a s » 0.5), but 
is consistent with the normalization, as ~ 0.9, implied by 
the cluster abundances. The normalization of our ACDM 
models, on the other hand, is in good agreement with both 
cluster and the 4-year COBE data. As explained above, 
the small difference in normalization between the ACDMi 
and ACDM2 was introduced to minimize differences be- 
tween the two power spectrum approximations used in 
these models. Normalization of the rCDM model is in- 
consistent with both the COBE normalization (as ~ 0.45) 
and with cluster abundance normalization (as ~ 0.52). 
This is motivated by our intent to use rCDM as a toy- 
model rather than a reasonable approximation to the real 
universe. The normalization of the rCDM is similar to 
normalizations of the CDM and ACDM models. There- 
fore, comparison between results of these models allows us 
to study effects of changing the shape of the power spec- 
trum and value of Qq. 

4. THE NUMERICAL SIMULATIONS 

4.1. Simulation parameters 

We have used two different A-body codes to carry 
out our simulations: the Adaptive Refinement Tree code 



(ART; Kra vtsov et al. 1997) and the AP 3 M code 7 ( pouch 
man 1991). Comparison of the results obtained with dif- 



ferent numerical codes allows us to insure the results are 
robust. AP 3 M code is an extension of the well-known 
P 3 M algorithm (Hockney & Eastwood 1981). The code 
performs hierarchical rectangular refinements in the high- 
density regions to reduce expensive particle-particle calcu- 
lations. The gravitational force is obtained by matching 
the gravitational forces calculated using FFT solver on 
the base and refinement grids and the small-scale force 
calculated using direct particle-particle summation (see 
Couchman 1991 for details). A total of four refinement 
levels were allowed during the course of the AP 3 M simu- 
lations presented here. The ART code also reaches high 
force resolution by refining all high-density regions with 
an automated refinement algorithm. The refinements are 
recursive: the refined regions can also be refined, each sub- 
sequent refinement having half of the previous level's cell 



7 The original public code was modified slightly to take into account Oq 7^ 1 models. 
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Table 2 
Parameters of simulations 



Code 


Model 


Run 


-^init 


^particle 

(h- l M Q ) 


-^stcps 


Resolution 
(kpc/h) 


Box 
(Mpc/h) 


-^part 


AP 3 M 


SCDM 


SCDM 


49 


3.5 x 10 9 


8000 


3.0 


30 


128 3 


AP 3 M 


OCDM 


OCDM 


109 


1.1 x 10 9 


7000 


4.7 


30 


128 3 


AP 3 M 


ACDMi 


ACDM 30 


64 


1.1 x 10 9 


4000 


3.0 


30 


128 3 


ART 


ACDM 2 


ACDM 60 


30 


1.1 x 10 9 


41300 


1.8 


GO 


256 3 


ART 


ACDM 2 


ACDM^RT 


45 


1.3 x 10 8 


13800 


0.9 


30 


256 3 


AP 3 M 


rCDM 


rCDM 


50 


3.5 x 10 9 


8000 


3.0 


30 


128 3 



size. This creates an hierarchy of refinement meshes of dif- 
ferent resolution covering regions of interest. The refine- 
ment is done cell-by-cell (individual cells can be refined 
or de-refined) and meshes are not constrained to have a 
rectangular (or any other) shape. This allows us to refine 
the required regions in an efficient manner. The crite- 
rion for refinement is local overdensity of particles: in the 
simulations presented in this paper the code refined an 
individual cell only if the density of particles (smoothed 
with the cloud-in-cell scheme; Hockney & Eastwood 1981) 
was higher than nth — 5 particles. Therefore, all regions 
with overdensity higher than S = n t h 2 3i /n, where fi is 
the average number density of particles in the cube, were 
refined to the refinement level L. For the two ART sim- 
ulations presented here, ACDM 60 and ACDM^ RT , n is 
1/8 and 1, respectively. The Poisson equation on the hi- 
erarchy of meshes is solved first on the base grid and then 
on the subsequent refinement levels. On each refinement 
level the code obtaines potential by solving a Dirichlct 
boundary problem with boundary conditions provided by 
the already existing solution at the previous level. There 
is no particle-particle summation in the ART code. The 
detailed description of the code is given in Kravtsov et al. 
(1997). Note, however, that the present version of the code 
uses multiple time steps on different refinement levels, as 
opposed to the constant time stepping in the original ver- 
sion of the code. The multiple time stepping scheme is 
described in some detail in Kravtsov et al. (1998; also see 
below). 

The information about the numerical parameters of the 
simulations is given in Table 2. The AP 3 M code was 
used to produce four simulations of different cosmologi- 
cal models 8 with the box size Lf, ox = 30/i _1 Mpc. The 
size of the box side is a compromise between requirements 
of the high spatial resolution (~ 2 — 4/i~ 1 kpc) and good 
statistics of halos. Nevertheless, for our most realistic 
model (ACDM) we also use the ART code to simulate 
a 60/i _1 Mpc box (ACDM 60 run). We estimate effects of 
the finite box size (see § 5.1) on our result s by compar 



ing results of lAOh *Mpc box simulation of Jenkins et al 
(1998 j with results of our 60/i _1 Mpc and 30/i 'Mpc boxes. 
All AP 3 M runs were done with 128 3 particles. The two 
ART runs used 256 3 particles. The initial conditions were 



set using the Zel'dovich approximation on uniform 128 3 , 
256 3 , and 512 3 meshes for AP 3 M runs, ACDM^ RT run, 
and ACDM 6 o run, respectively. The seed used to gener- 
ate the Gaussian random density field was the same in all 
of our AP 3 M runs, but different for each of the two ART 
runs. All of the simulations are started at the moment 
of time when the rms density fluctuations at the Nyquist 
wavelength \N yq are still linear: o{\Nyq, Zi) ~ 0.1 — 0.2. 
All AP 3 M runs evolved during a period in which the linear 
growth factor increased by a factor of 50. This explains 
different values of Zi„ lt for different models in Table 2. 

As was explained in §1, the purpose of our study was 
to compute the correlation function and the bias account- 
ing for all DM halos, including those inside groups and 
clusters. To assure that halos do survive in clusters the 
force resolution should be ~ 1 — 3/i _1 kpc (Moore et al. 
1996; KGKK). Furthermore, if we aim to study galaxy- 
size halos, the mass resolution should be <; 10 9 Ii~ 1 Mq to 
resolve galaxy-size halos (M J> 10 11 /i _1 Mq) with at least 
ss 100 particles. The compromise between these considera- 
tions and the computational expense determined the force 
and mass resolution of our simulations (see Table 2). The 
ART code integrates the equations of motion in comoving 
coordinates. However, the refinement strategy of the ART 
code is designed to effectively preserve the initial phys- 
ical resolution of the simulation (see below). The peak 
resolution is reached by creating a refinement hierarchy 
with six levels of refinement. In the AP 3 M runs the force 
resolution n (spline softening length) was kept constant in 
comoving coordinates while fluctuations are still in the lin- 
ear regime and is then set to be constant in physical units. 
The switch occurs at the moment when first galaxy-size 
halos start to collapse (z ~ 5 — 10) for our simulations. 
We chose to maintain fixed comoving resolution until it 
reaches ~ 3/i -1 kpc (physical) (at ~ 5 — 10). At later mo- 
ments the resolution is fixed to this value in physical co- 
ordinates (the exception is the OCDM model in which the 
resolution was set to 4.7ft.~ 1 kpc by mistake). The dynam- 
ical range Lt, ox /r] of the simulations implied by the force 
resolution is « 16,000 (32,000 formal) for the ART runs 
and 6000 - 10000 for the AP 3 M runs. The dynamic range 
of the AP 3 M runs is just enough to keep the initial physi- 
cal resolution («2- 5/i _1 kpc). The ART code integrates 



'Runs SCDM, ACDM30 , and tCDM have not been completed due to the high computational expense. The simulations were stopped at 
: 0.3. 
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following we describe specific parameters of the BDM used 
to construct halo catalogs used in our study (the main pa- 
rameters are listed in Table 3). 

The radius of a halo assigned to it by the algorithm is ei- 
ther its virial radius 10 or lSO/i^kpc, whichever is smaller. 
The latter is approximately the maximum virial radius we 
would expect for a galaxy-size halo. The mass and ra- 
dius are very poorly defined for the satellite halos due to 
the tidal stripping which alters a halo's mass and physical 
extent. Therefore, in this study we will use maximum cir- 
cular velocity V max as a proxy for halo mass. This allows 
us to avoid complications related to the mass and radius 
determination for satellite halos. For isolated halos V max 
and the halo's virial mass are directly related. For exam- 
ple, for a halo with a density distribution described by the 
Navarro, Frenk & White (hereafter NFW; 1996) profile 
p(r) (x x _1 (l + x)~ 2 (x = r/R s ; R s is scale-radius): 



V. 



2 

max 



GAL, 



/(2) 



Rvir /(c) 2 



(9) 



the evolution in comoving coordinates. Therefore, in order 
to prevent degradation of force resolution in physical coor- 
dinates, the dynamic range between the start and the end 
(z = 0) of the simulation should increase by (1 + 2,): i.e., 
for our simulations reach 512 x (1 + zi) = 15,872. This 
is accomplished with the prompt successive refinements 
during the simulations. 

The time stepping of the AP 3 M and ART codes is rather 
different. First of all, the codes integrate the equation 
of motion using different time variables: the time in the 
AP 3 M code and the expansion factor in the ART code. In 
the AP 3 M runs the time step is constant and is the same 
for all particles. In the ART runs, as was noted above, the 
particles residing on different refinement levels move with 
different time steps. The particles on the same level, how- 
ever, move with the same step. The refinement of the time 
integration mimicks spatial refinement and step for each 
subsequent refinement level is two times smaller than the 
step on the previous level. The global time step hierarchy 
is thus set by the step Aao at the zeroth level (uniform 
base grid). The step on level L is then Aa^ = Aeto/2 L . 
The choice of an appropriate time step for a simulation 
is dictated by the adopted force resolution. The number 
of time steps in our simulations is such that the rms dis- 
placement of particles during a single time-step is always 
less than rj/4 (less than 1/4 of a cell in the ART code) 9 . 
The size of the time-step, At, for the AP 3 M runs was cho- 
sen to be sufficiently small to sati sfy the stability criteria 
of the numerical integration (e.g., Efstathiou et al. 1985) 
throughout the entire run. In the case of ART runs, the 
value of Aao — 0.0015 was determined in a convergence 
study using a set of smaller 64 3 particle simulations de- 
scribed in Kravtsov et al. (1998). In both AP 3 M and ART 
runs the energy was conserved with an accuracy <; 1%. 

4.2. Identification of halos 

Identification of DM halos in the very high-density en- 
vironments (e.g., inside groups and clusters) is a challeng- 
ing problem. Traditional halo finding algorithms, such as 
fricnds-of- friends (e.g., Davis et al. 1985) or "overdensity- 
200" (e.g., Lacey & Cole 1994), cannot be used. These al- 
gorithms are not designed to search for substructure; they 
identify an isolated halo above virial overdensity as a single 
object and cannot account for the internal substructure. 
Our goal, however, is to identify both isolated halos and 
halos orbiting within larger systems ( "sub- halos" ) . The 
problems associated with halo identification within high- 
density regions are discussed in KGKK. In this study we 
use a halo finding algorithm called Bound Density Maxima 
(BDM; see KGKK). A detailed description of the working 
version of the BDM algorithm used here can be found in 
Klypin & Holtzman (1997). Other recently developed al- 
gorithms capable of identifying satellite halos are described 
in Ghigna et al. (1998) and KGKK. The main idea of the 
BDM algorithm is to find positions of local maxima in 
the density field smoothed at a certain scale and to apply 
physically motivated criteria to test whether the identified 
site corresponds to a gravitationally bound halo. In the 

9 Note, however, that the distance traveled by the fastest moving particle in one time-step in AP 3 M runs can be larger than r;, especially at 
late times and in the ACDM30 run. In the ART code particles do not move further than ~ 0.5 cells in a single time step, where the cell size 
and time step for particles located on the refinement level L are Axo/2 L and Aao/2 L , respectively. 

10 The virial overdenisty 5th is set in accord with the prediction of the top-hat collapse model. Note that &th depends on the cosmological 
model (e.g., Kitayama &; Suto 1996). 



where M v i r and R V i r are virial mass and radius, fix) = 
ln(l+x) — x/(l+x), c = R vir /R s . While for the sub- halos 
V m ax may not be related to mass in any obvious way, it 
is still the most physically and observationally motivated 
halo quantity. The limiting radius of 150/i _1 kpc is suf- 
ficient to determine the V max for galaxy-size halos. The 
cluster-size halos are not explicitly excluded from the halo 
catalogs. We assume therefore that the center of each clus- 
ter can be associated with a central cluster galaxy. The 
latter (due to the lack of hydrodynamics and other rele- 
vant processes) cannot be identified in our simulations in 
any other way. 

The density maxima are identified using a top-hat fil- 
ter with radius r s ("search radius"). The search is per- 
formed starting from a large number of randomly placed 
positions ("seeds") and proceeds by moving the center of 
mass within a sphere of radius r s iteratively until conver- 
gence. In order to make sure that we use a sufficiently 
large number of seeds, we used the position of every tenth 
particle as a seed. Therefore, the number of seeds by far 
exceeds the number of expected halos. The search radius 
r s also defines the minimum allowed distance between two 
halos. If the distance between centers of any of the two 
halos is < 2r s , only one halo (the more massive of the two) 
is left in the catalog. A typical value for the search radius 
is (10— 20)/i _1 kpc. We set a lower limit for the mass inside 
the search radius M{< r s ): halos with M(< r s ) < M m ;„ 
are not included in the catalog. This is done to exclude 
pure poisson fluctuations from the list of halo candidates. 
Some halos may have significant substructure in their cores 
due, for example, to an incomplete merger. Such cases ap- 
pear in the catalogs as multiple (2-3) halos with very simi- 
lar properties (mass, velocity, radius) at small separations. 
Our strategy is to count these as a single halo. Specific cri- 
teria used to identify such cases are: (1) distance between 
halo centers is <; 150/i~ 1 kpc, (2) their relative velocity in 
units of the rms velocity of particles in the halos Ad jv is 
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Table 3 

Parameters of Halo Finding Algorithm 



Run 


5th 


M min 


T 








(kpc/h) 


SCDM 


180 


7 x 10 10 


13 


OCDM 


180-400 


2 x 10 10 


13 


ACDM 30 


180-340 


2 x 10 10 


13 


ACDM 60 


180-340 


10 10 


20 


ACDM^fT 


180-340 


10 9 


10 


rCDM 


180 


7 x 10 10 


13 



less than 0.15, and (3) the difference in mass is less than 
factor 1.5. Only the most massive halo is kept in the cat- 
alog. 

It is obvious that for a statistical study it is important 
to be confident that our BDM algorithm does not miss a 
large fraction of halos. While a substantial effort is made 
to reject fake halos, it is important to make sure that all 
real halos are included in the catalog. The algorithm is 
very efficient at finding isolated halos and the major dif- 
ficulties are in identification of halos in crowded regions. 
Halo identification in such regions is complicated by the 
large number of halos and by the high-density background 
of fast moving particles. Therefore, performance of the 
algorithm in such regions is a good indicator of its overall 
performance. In fact, the parameters of the halo finder 
used to construct our halo catalogs were tuned by visual 
inspections of the most difficult and complicated regions. 
An example of such regions is shown in Figure 1 . The fig- 
ure shows the distribution of the DM particles in a Virgo- 
type cluster 11 in the ACDM model (ACDM 60 run). To 
enhance the substructure, particles are color-coded on a 
grey-scale with a rendering algorithm based on the local 
particle density (see caption). A large number of distinct 
and compact DM halos is clearly present within the virial 
radius of the cluster. Figure 2 shows positions of DM halos 
identified by the BDM code in the same volume. There 
are 121 halos in the plot, each halo having more than 20 
bound particles. More stringent limits of more than 30 
particles and a maximum circular velocity, V ma x, larger 
than 80km/s produce 98 halos. All distinct halos visible 
in Figure 1, are identified. 

4.3. Survival of halos in clusters 

In § 1 we have stressed that the main new feature 
of the analysis presented in this paper is the identifica- 
tion and use of halos located inside the virial radius of 
other halos (satellites or sub- halos). Particularly, one of 
the main goals is to include halos inside cluster-size ha- 
los. This would allow an estimation of the halo-halo cor- 
relation function down to unprecedentedly small scales 
(~ 150/i _1 kpc). However, we can consider this estimate 
robust only if we are confident that no halos (or a very 

11 The distribution of matter and halos in this cluster is discussed in 
12 similar to those of KGKK but with all parameters appropriate for 
13 Here we quote all resolutions for the spline kernel rj of the AP 3 M, 



small fraction of them) arc artificially destroyed in clusters 
due to the insufficient resolution. KGKK have discussed 
analytical estimates and numerical experiments that could 
be used to address this issue. Following the approach 
of KGKK, we have run a series of small TV-body simu- 
lations 12 using the direct-summation Aarseth's code (Bin- 
ney & Tremain 1987). The basic setup of the simulations 
is as follows. A DM halo of virial mass lO 12 ft _1 M (con- 
taining a few thousand particles within the virial radius) is 
constructed. The initial equilibrium density profile of the 
halo is described by the NFW formula. The halo is then 
placed on an orbit in a constant potential corresponding to 
a galaxy cluster of mass 2 x lO 14 ft~ 1 M0with the NFW den- 
sity distribution. The particular numbers quoted here are 
intended to mimic the orbital evolution in the cluster pre- 
sented in Figs. 1-3. The orbital evolution of the halo was 
studied for different orbits and different (mass and force) 
resolutions. Both the mass and the force resolution were 
varied by more than a factor of ten. The orbital evolution 
of mass bound to the halo converges at the force resolu- 
tion 13 m 3 — 4/i _1 kpc: i.e., the evolution of bound mass in 
runs with higher resolution is identical. In fact, during the 
first i=a 5 Gyrs the resolution of w lO/i kpc is adequate. 
The mass loss in this case is somewhat higher which leads 
to total destruction after 5 Gyrs, whereas resolution of 
3/i _1 kpc allows halo to survive during the Hubble time. 
These experiments have also shown that mass resolution 
of w 10 9 /i _1 M Q is sufficient, provided that force resolu- 
tion is high. Two runs with force resolution of 3/i _1 kpc 
and with mass resolutions of 10 8 /i~ 1 M Q and 10 9 /i~ 1 M Q 
resulted in identical mass evolution. It is worth noting 
that in these experiments the halo was followed until it 
was totally destroyed by the tidal field (;> 5 Gyrs in most 
cases). In real simulations, however, clusters form only at 
z <; 1, and most of the accreted halos spend <J 5 Gyrs in 
clusters. During this time halos lose 80-90% (depending on 
the orbit) of their initial mass. Thus, if the mass resolution 
is w 10 9 /i _1 M Q , the halos with initial mass ;> 10 11 h~ 1 M Q 
can be identified even after spending a substantial time in 
a cluster. 

Details of the evolution of halos in clusters depend sen- 
sitively on the parameters of the halo orbit. Specifically, 

the next section, 
the ACDM model. 

the corresponding resolution is rj/2 for the ART. 
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Fig. 1. — Distribution of dark matter in a Virgo- like clus- 
ter in the ACDMeo simulations. The cluster virial mass is 
2.45 x 1O 14 /i _1 M0 and the corresponding 3D velocity dispersion 
is « 1022km/s. The particles inside a sphere of the radius of 
1.5ft _1 Mpc (solid circle) are shown. The size of the small box, 
shown to provide the comparison scale, is 100/i _1 kpc. To en- 
hance the contrast, we have color-coded DM particles on a grey 
scale according to their local density: intensity of each parti- 
cle is scaled as the logarithm of the density difference pis — P75, 
where the densities were obtained using top-hat filter with radii 
15/i _1 kpc and 75/i -1 kpc. 

the mass loss rate depends on the pericenter and on the 
eccentricity of the orbit. A halo on a very eccentric orbit 
survives for a considerably longer time than a halo on a 
circular orbit, even if the radius of the latter is larger than 
the pericenter. For example, a halo on a circular orbit 
with a radius of 250/i~ 1 kpc was totally destroyed in less 
than 5 Gyrs, while a halo on a very eccentric orbit with 
the pericenter of only 125/i _1 kpc survived for more than 
10 Gyrs. The explanation for this is simple: a halo on a 
circular orbit spends all of its time in a high-density clus- 
ter core suffering a steady mass loss, while a halo on the 
eccentric orbit spends only a small fraction of its orbital 
period in the core. 

With the resolution quoted above, only a relatively small 
fraction of halos can be tidally destroyed in clusters. Only 
the halos with the apocenter <^ 300/i _1 kpc are subject to 
the destruction. It can be expected that halos which are 
accreted when the cluster was young and its radius was 
small (~ 300 — 500/i _1 kpc) have small apocenters. The 
fraction of such halos in a z — cluster is small. Our es- 
timates show that for halos with large apocenters the dy- 
namical friction cannot bring the apocenter considerably 
closer to the cluster center. The tidal stripping reduces 




Fig. 2. — The region of the space shown in Figure 1 with cir- 
cles representing the DM halos found by the BDM halo finder. 
Area of each circle representing a halo is proportional to the 
halo's maximum circular velocity. There are 121 halos in the 
plot, each of which has mass > 2 x lO lo /i _1 M0 (more than 20 
bound particles). 



the halo mass very efficiently, thus increasing the friction 
time. 

The tidal stripping that halos suffer in clusters has an 
important effect on their density profiles. It appears that 
the halo profile is affected at all radii and not only at radii 
close to the tidal radius. Indeed, the trajectories of parti- 
cles in a halo are mostly eccentric and after a crossing time 
the absence of the stripped particles will be "felt" by the 
whole halo. The change in the density profile means, of 
course, a change in the maximum circular velocity of the 
halo. The effect is not dramatic because the fraction of 
particles on highly eccentric orbits is rather small and the 
central regions of halos are thus affected the least. Typ- 
ically, the circular velocity of the halos is reduced during 
their evolution in the cluster by 20%-30%. This correc- 
tion was taken into account in Figure 3 by reducing the 
velocity cut from 120 km/sec at z=1.0 to 100 km/sec at 
z=0.0. 

The tidal destruction is important only for cluster-size 
halos (where tidal fields are the strongest). There are 
three massive clusters in our ACDM 6 o simulation with 
mass ^ 2 x lO 14 /i _1 M . There is one "Coma clus- 
ter" with velocity dispersion a^o = 1654 km/s and mass 



10 



COLIN ET AL. 



M = 6.4 x 1O 14 /i _1 M inside 1.43/i _1 Mpc radius. The 
cluster contains 201 halos of mass > 3 x 10 10 /i -1 Mq inside 
a 2/i _1 Mpc box and 9 halos inside 0.3/i _1 Mpc. Unfortu- 
nately, the cluster has suffered a recent major merger and 
it was difficult to study the radial distribution of halos. 
One of the other two clusters also shows an indication of 
an ongoing merger. It has two sub-clusters in the cen- 
tral region separated by w 0.5/i _1 Mpc. Therefore, we will 
focus on the third cluster which has a relatively regular 
appearance (see Fig.l). 

The cluster contains 121 bound halos with > 20 particles 
inside 1.5fo — 1 Mpc. There are 231,200 dark matter parti- 
cles inside the virial radius of 1.28/i _1 Mpc. The virial mass 
of the cluster has increased from M v ; r = 7.9 x lO 13 ft. _1 M0 
at z = 1 to M vir = 2.4 x lO 14 /i _:l M at z = 0. Figure 
3 shows the density profile of matter and the radial num- 
ber density profile of DM halos inside the cluster. The 
DM density profile at z = is well approximated by the 
NFW profile. For comparison, we also present the density 
profile of this cluster at redshift one. Note that at z = 1 
the radius is given in proper units and the mean matter 
density (po) is estimated at z = 0. Similarly, the bottom 
panel of Figure 3 shows the number density profile of ha- 
los in the cluster at z — and at z = 1. Halos with more 
than 30 bound particles and with limits on the maximum 
circular velocity V max > lOOkm/s and V max > 120km/s 
(the change is for the reasons explained above) for z = 
and z = 1, respectively, were used. The mean number 
density of halos, no, was estimated using all halos in the 
simulation within the above velocity limits (N^ = 7628 at 
z = 0). The profile at z = 1 is also rescaled into proper 
units. Figure 3 clearly indicates that the number of ha- 
los in the central 300/i -1 kpc (proper radius) has declined 
substantially from z = 1 to the present epoch. In the cen- 
tral 300/i _1 kpc there are three times as many halos (24) 
at z = 1.0 than at z = 0. If we interpret the difference as 
due to tidal destruction, then 16 halos were destroyed in 
the central part. The situation is different at larger radii: 
there are 57 and 50 halos at 0.3 < r < 1.3/i~ 1 Mpc for 
z = 1 and z = 0, respectively. Note, that we compare the 
number of halos in the same proper volume at these two 
moments: the volume corresponding to the virial radius of 
the cluster at z = (the virial radius at z = 1 is smaller). 
The number density profile is thus virtually unaffected at 
r £ 300/1-^0. 
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Fig. 3. — The density profile of the cluster shown in Fig- 
ures 1 and 2. Top panel: Dark matter density in units of 
the mean matter density at z = (solid line) and at z = 1 
(dot-dashed line); the dashed line shows the Navarro, Frenk 
& White (1996) fit to the z = density profile. The z = 1 
density profile is given in proper units: the radius is given 
in proper scale and the mean matter density is estimated at 
z — 0. Bottom panel: Number density profiles of halos in the 
cluster at z — (solid circles) and at z — 1 (open circles) as 
compared to the z — density profile of dark matter (solid 
curve). The error-bars show la poisson errors. Halos with 
more than 30 bound particles and with maximum circular ve- 
locity larger than lOOkm/s at z — and larger than 120km/s 
at z = 1 were used to estimate the average density no (see 
§ 4.3 for details). The profile at z = 1 is rescaled into proper 
units similarly to that of the dark matter in the top panel. 



5. RESULTS 

In this section we will present results on the evolution 
of the 2-point correlation function and bias (as defined 
by eq. [1]) in the simulations described in the previous 
section. The presentation of the results is split in two 
subsections. In the first of them we present the results of 
our largest simulation: ACDMgo- In the second section 
we present the results of the rest of our simulations of 
different cosmological models. Thus, we will first focus 
on the results of the most realistic of the studied models, 
ACDM (see § 3), and then will discuss the differences be- 
tween cosmological models. However, before we proceed 
with the presentation of the results, we will first compare 
our estimate of the dark matter correlation function with 
estimates which have been performed by different authors 



and with different numerical codes. 

5.1. The dark matter correlation function in the ACDM 
model: comparison with other studies 

In Figure 4 we compare the 2-point correlation function 
of the dark matter in the ACDM model (the ART run 
ACDMgo) with similar estimates presented by Klypin et 
al. (1996) and by Jenkins et al. (1998). The two latter 
estimates were done with different codes (PM in Klypin 
et al.; and AP 3 M in Jenkins et al.) and different reso- 
lutions (cell size of 62/i _1 kpc for the PM, and 30/i _1 kpc 
(Plummer) for the AP 3 M). All simulations followed 256 3 
particles, although physical mass resolution was different 
due to the different box sizes (50/i _1 Mpc, 60/i~ 1 Mpc, and 
141.3/i -1 Mpc for the PM, ART, and AP 3 M simulations). 
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Fig. 4. — The comparison of the correlation functions of 
the dark matter £dm in the ACDM model estimated by dif- 
ferent authors with different numerical resolutions and codes. 
The solid curve shows in our ACDMgo run, simulated us- 
ing the ART code. The dashed curve shows £<jm estimated 
by Jenkins et al. (1998) using the AP 3 M code. The crosses 
show £dm estimated by Klypin et al. (1996) using the PM 
code. The AP 3 M and PM correlation functions were rescaled 
as described in § 5.1 to account for the difference in normal- 
ization between different simulations. All simulations followed 
evolution of 256 3 particles. The vertical lines indicate formal 
force resolution for each code (the line for the ART code at 
1.8h~ kpc is off the plot). 

Although the cosmological model was exactly the same in 
all three estimates 14 , the normalization of the power spec- 
trum was slightly different. The rms mass fluctuations on 
scale of 8/i- x Mpc a 8 was 0.9, 1.0, and 1.1 for the AP 3 M, 
ART, and PM simulations, respectively. Therefore, we 
have multiplied (divided) the DM correlation function of 
AP 3 M (PM) simulation by l.l 2 in order to account for the 
differences in normalization. Figure 4 indicates that there 
is a very good agreement between all estimates at scales 
« (0.2 - 2)/i" 1 Mpc. The ART and AP 3 M estimates agree 
to better than 10% at scales 0.03 — 7/i _1 Mpc! On larger 
scales the ART correlation function has a smaller ampli- 
tude than that of the AP 3 M due to a factor of 2.35 smaller 
box size of the former. This result is in agreement with 
previous studies (e.g., Colin et al. 1997 and references 
therein) which concluded that the correlation function is 
underestimated at scales 0.1 of the simulation box size. 
Nevertheless, the agreement is striking at smaller scales, 
given all the differences (including cosmic variance) be- 
tween the estimates. We conclude, therefore, that the 
correlation functions presented in the next subsection are 
robust at scales <; 7/i~ 1 Mpc. This scale is, of course, lower 
(« 3 — 4/i _1 Mpc) for the 30/i _1 Mpc simulations presented 
in § 5.3. 

14 Note, however, that there are small differences in the approxima 
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Fig. 5. — The evolution of the 2-point correlation function 
of the dark matter (bottom panel) and halos (top panel) in 
the ACDM60 simulation. Only halos with maximum circular 
velocity > 120km/s were used to estimate the halo correlation 
functions. Poisson errors for the halo correlation functions are 
negligible at scales :> 0.2/i -1 Mpc and are not shown for clar- 
ity; at scales < 0.2/i _1 Mpc the errorbars are <J 20% (see § 5.2 
for details). The best fit parameters of the power- law fit to 
the halo correlation functions are given in Table 4. 

5.2. Evolution of the correlation function and bias in the 
ACDM model 

In Figure 5 we plot the evolution of the correlation func- 
tion of both the dark matter ^ TO and the DM halos £,hh 
in our ACDMgo simulation. All results in this and the 
following sections are presented in comoving coordinates. 
The halo-halo correlation function was constructed using 
halos with V ma x > 120km/s. Note that although the ear- 
liest moment at which we show the correlation function is 
z = 5, the first halos in this simulation collapsed at redshift 
z f« 10 and hundreds of halos are identified at z k, 7—8. 
The number of halos in the V max > 120km/s catalog is 
approximately 4300, 10, 000, and 7500 for redshifts of 5, 2, 
and 0, respectively. The good statistics result in a very ac- 
curate estimate of the correlation function. The number of 
pairs per scale bin used to estimate ^hh is > 50 in all cases. 
Typically, 60 - 150 for r < Q.2h~ 1 Mpc, 200 - 1000 for 

used for the power spectrum between the simulations. 
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Table 4 

Power-law fits to £hh IN THE ACDM 60 simulation 



V m ax > 120km/s 



V m ax > 150km/s 



Vmax > 200km/s 



7 



ro 



7 



ro 



7 



0.0 
0.5 
1.0 
1.5 
2.0 
3.0 
4.0 
5.0 



4.864±0.011 
4.145±0.009 
3.760±0.007 
3.419±0.006 
3.103±0.006 
2.974±0.006 
3.206±0.009 
3.602±0.013 



1.704±0.004 
1.684±0.004 
1.694±0.003 
1.642±0.003 
1.588±0.003 
1.527±0.004 
1.516±0.005 
1.534±0.007 



4.789±0.024 
4.040±0.019 
3.944±0.015 
3.823±0.013 
3.831±0.011 
3.453±0.011 
3.356±0.011 
3.661±0.016 



1.687±0.009 
1.726±0.009 
1.762±0.007 
1.715±0.006 
1.629±0.006 
1.536±0.006 
1.498±0.006 
1.516±0.008 



5.082±0.062 
4.534±0.049 
4.387±0.038 
4.288±0.035 
4.375±0.036 
4.327±0.037 
4.733±0.044 
4.829±0.053 



1.650±0.022 
1.777±0.020 
1.869±0.016 
1.831±0.016 
1.702±0.016 
1.683±0.016 
1.598±0.017 
1.578±0.020 



a Column description: z is redshift; r and 7 are the best fit parameters of the power-law 
£,hh{r) — (r/r )~ 7 fit to the comoving halo-halo correlation function; r is comoving in units /i _1 Mpc. 



0.2 < r < l/i^Mpc, and > 1000 for the larger scales. The 
pure poisson errors associated with each of the points are 
thus neglegibly small, except for the first 2-3 bins (where 
poisson errors are still small: <J 20%). There are subtler 
errors associated with radial binning, but these are also 
less than a few per cent. 

Figure 5 shows that the shapes of the matter and 
halo-halo correlation functions are quite different. The 
matter correlation function changes its shape from al- 
most a power-law to a complicated shape. The slope of 
£dm at scales £ 0.5/i _1 Mpc stays approximately constant 
throughout the evolution, while at the larger scales ^ m 
significantly steepens. The amplitude of the DM correla- 
tion function increases from z = 5 to z = by factors of 
s=y 60 and w 10 at small and large scales, respectively. The 
halo-halo correlation function behaves very differently Its 
shape can be well-described by a power-law at all epochs. 
The amplitude of ^hh evolves non-monotonically: it de- 
creases somewhat from z = 5 to z — 3, and then gradually 
increases. The evolution, however, is much more modest 
than the dramatic evolution of the £d m : the maximum dif- 
ference in amplitude among any two epochs is only a factor 
of two. The details of the £hh evolution are illustrated in 
Figure 6. This figure shows the evolution of the £,hh ampli- 
tude at a variety of different comoving scales (indicated on 
the right) for the three ACDM runs from our set of simu- 
lations (see Table 2). Different initial conditions of these 
runs allow us to evaluate the cosmic variance, while differ- 
ent particle masses and spatial resolutions (by a factor of 8 
and 2, respectively) of ACDM 60 and ACDMg RT allow us 
to check for the resolution effects. Comparison of £,hh{r, z) 
in the 30ft -1 Mpc runs simulated with different codes and 
with significantly different resolutions shows that the two 
runs agree very well. The KCT)M^ T simulation had the 
best mass and force resolution, yet, we do not find any 



visible systematic differences at scales <J 2ft~ 1 Mpc with 
the other two simulations. There is an indication that at 
r ^ 2/i~ 1 Mpc the amplitude in the 30/i -1 Mpc simulations 
is systematically lower than that in the ACDMgo, which 
can be explained by the effects of finite box size. This is 
in agreement with expectation that the finite size effects 
become important at scales > 3/i _1 Mpc (w 0.1 of the box 
size). At smaller scales, where we expect the SO/i^Mpc 
simulations to produce correct results, the agreement is 
very good 15 . Besides the illustration of a very mild evo- 
lution of £hh, Figure 6 also shows that there is a common 
feature of the evolution. Although the exact evolution 
depends to some extent on the scale, the amplitudes at 
all scales are quite high (as high or higher as they are 
at z = 0) at very high redshift (z s=a 7). The ampli- 
tude then decreases until z ~ 2 — 4, and grows steadily 
at lower redshifts. It is important to note that this evo- 
lution is more complicated than simple evolution models 
often used in the observational and theoretical analyses: 
£(r,z) = (r/ro)~ 7 (l + z)-( 3 ~ 7+e ). Figure 6 shows that 
parameter e estimated by such analyses would depend not 
only on the redshift range used, but also on the scale at 
which the amplitude is measured (as well as on other pa- 
rameters such as the object's mass). This calls into ques- 
tion the usefulness of such a simplistic approach (see also 
arguments in Moscardini et al. 1998). Note that there is 
also some observational evidence (Giavalisco et al. 1998) 
indicating that the above parameterization is a poor de- 
scription of the observed galaxy clustering evolution. For a 
limited range of redshifts z ~ — 1, we find only very weak 
evolution of halo clustering in comoving coordinates indi- 
cating a value of e w — 1. This value seems to be favored 
by observations of galaxy clustering at these redshifts (e.g., 
Postman et al. 1998). 



5 Note that some cosmic variance is expected for these box sizes. 
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Fig. 6. — The evolution of the correlation function of halos 
in ACDM60 simulation at various comoving scales. Each curve 
shows the amplitude of the correlation function at a fixed co- 
moving radius indicated on the right side of the panels. The 
correlation functions were estimated for halos with circular 
velocity larger than 150km/s (top) and > lOOkm/s (bottom). 
The solid curves and open circles indicate amplitudes in the 
ART 60ft _1 Mpc and 30/i~ 1 Mpc simulations, respectively; the 
triangles show the amplitude in the 30/i _1 Mpc AP 3 M simula- 
tion (see § 5.2 for details and discussion). 

The power- law fits to £hh of the form £,hh(r) = (r/ro)~ 7 
for various epochs and for halo catalogs with different cuts 
in the maximum circular velocity are presented in Table 
4. Direct (rather than linear fits to \g£hh ~ lg r ) weighted 
power-law fits were done to £hh with the Levenberg- 
Marquardt method described in Press et al. (1992). Each 
bin of the correlation function was weighted with its Pois- 
son error <j£ — (^hh + l)/-\7^p> where n p is the number of 
halo pairs in the bin. Visual inspection shows that power- 
law fits are very successful at scales ;> 0.3/i~ 1 Mpc, while 
at smaller scales there are « 20 — 30% deviations in some 
cases. The goodness of the fits is represented in rather 
small formal errors of the best fit parameters r and 7. We 
have estimated how these parameters change if the corre- 
lation function is re-binned differently and found that the 
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Fig. 7. — Bottom panel: Comparison of the halo correlation 
function in the ACDM60 simulation with the correlation func- 
tion of the APM galaxies (Baugh 1996) . Results for halos with 
maximum circular velocity larger than 120km/s, 150km/s, and 
200km/s are presented by the solid, dot-dashed, and dashed 
curves, respectively. The dotted curve shows the dark matter 
correlation function. Note that at scales > 0.3/i -1 Mpc the 
halo correlation function does not depend on the limit in the 
maximum circular velocity (see § 5.2 for details). Top panel: 
Dependence of bias on scale and maximum circular velocity. 
The curve labeling is the same as in the bottom panel, ex- 
cept that the dotted line now represents the bias of halos with 
Vmax > 100km/s. 

change is always <J 3%. The examination of the Table 4 
shows that typical ranges of ro and 7 are ~ 3 — 5/i _1 Mpc 
and « 1.5 — 1.7, respectively. For all halo catalogs, param- 
eters evolve slowly with redshift. The correlation length 
ro decreases somewhat between redshifts of 5 and 3, and 
then increases steadily until z = 0. The evolution of 
7 is even slower with the tendency for 7 to increase by 
~ 10% from z = 5 to the present epoch. An important 
and interesting point is that the correlation amplitude, 
and hence the value of ro, are quite different for halo cat- 
alogs with different V max cuts. The correlation lengths 
for Vmax > 120km/s and V max > 200km/s catalogs differ 
by 25% at z > 3, while the difference is only « 4% at 
z = 0. This means there is a mass segregation of halo 
clustering properties at high redshift, which, however, is 
erased during the subsequent evolution. 

During the last few years, there has been tremen- 
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dous progress in the observational studies of high-redshift 
galaxy clustering. We will discuss how our results on the 
halo clustering evolution compare with the results of ob- 
servations in §6. Here, however, we present a comparison 
of the £hh with the most accurate measurement of the 
galaxy correlation function £ gg (at z w 0) made using the 
APM galaxy survey (Baugh 1996). Figure 7 shows the 
z = correlation functions of halos and dark matter in 
the ACDM model and the real space APM galaxy cor- 
relation function. The halo-halo correlation function was 
estimated for halo catalogs with cuts in the maximum cir- 
cular velocity of 120km/s, 150km/s, and 200km/s. The fig- 
ure shows striking agreement between the halo and galaxy 
correlation functions: at scales ;> 0.3h~ 1 Mpc the corre- 
lation functions of all halo catalogs match both the shape 
and the amplitude of the £ gs . The correlation function 
for V max > 150km/s catalog agrees with APM £ 9S within 
errors at all probed scales. As we noted above, the dif- 
ferences that exist between the catalogs at high redshifts 
virtually vanish during the course of the evolution. This 
is manifested in the similarity of £hh for different catalogs 
on scales ;> 0.3/i _1 Mpc. Note, however, that this does 
not mean that ^hh has no mass dependency. Rather, the 
result means that by z = any mass dependence of the 
correlation function vanishes when averaged over a range 
of galactic masses. As was explained above, the poisson er- 
rors of the halo correlation functions shown in Figure 7 are 
very small and were not shown for clarity. The robustness 
of the result can be estimated, however, by comparing 
of Vmax > 120km/s and V max > 150km/s catalogs. The 
number of halos in these catalogs is significantly different: 
4708 and 2480, respectively. This makes the halo sam- 
ples largely independent. The correlation functions agree, 
however, within the poisson errors. 

While the correlation function of halos matches that of 
galaxies very accurately, the correlation function of mat- 
ter {.dm matches £ gg neither in shape nor in amplitude 16 . 
The amplitude is matched only at scales ;> 4 — 5/i Mpc. 
At smaller scales it is much higher than the amplitude of 
the APM £ gg , implying that DM halos are anti-biased at 
these scales with respect to the dark matter. Moreover, 
the difference in shape between ^hh and ^ m implies that 
the bias is scale-dependent. The scale dependence of the 
bias {^J S,hh{f) / £,dm{f)) for the halo correlation functions 
is shown in the top panel of figure 7. The bias varies 
significantly at scales <; 5ft, Mpc in the range ~ 0.5 — 1. 
Moreover, as was shown in Figure 5, the shape of the 
correlation function of dark matter differs from that of 
£hh at all epochs and evolves much more strongly than the 
correlation function of halos. The former fact implies that 
the bias is scale- dependent at all epochs, while the latter 
means that the bias evolves rapidly with cosmic time. The 
evolution of bias in the ACDMgo run is illustrated in Fig- 
ure 8. The evolution is shown for different halo catalogs 
and at different scales. The bias evolves very rapidly from 
value of ~ 3 — 5 at z ss 5 to ~ 0.5 — 1 at z = 0. The evolu- 
tion depends on the velocity (or mass) cut of the catalog 
at high z: the halos in the catalogs with higher velocity 
cuts exhibit stronger clustering. This difference vanishes, 
however, at z <; 0.5. The evolution of the scale dependence 

le This result is in agreement with conclusions of Jenkins et al. (1998). 
that calculated by Jenkins et al. 




0.4 1 1 1 1 1 L 

1 2 4 6 



(1+z) 

Fig. 8. — Top panel: The evolution of bias at comoving 
scale of 0.54ft" 1 Mpc for halos with different lower limit on the 
maximum circular velocity in the ACDMeo simulation. The 
bottom panel shows dependence of the bias on {comoving) scale 
for halos with maximum circular velocity > lOOkm/s. 

of the bias is also interesting. At high redshifts the bias 
was larger at small scales. At small redshifts the halos are 
almost unbiased (b 1) on a few megaparsec scales and 
are anti-biased (b « 0.5 — 0.6) on small scales. 

5.3. Evolution of the correlation function and bias in 
different cosmological models 

Is evolution of bias observed in the ACDM simulations 
specific to this model or is this evolution similar for all 
of the models? We address this question by compar- 
ing results presented in the previous section with results 
of the 30/i _1 Mpc simulations of other cosmological mod- 
els (see Table 2). Figures 9—12 show the correlation 
functions of halos ^hh and the dark matter ^ m for the 
four AP 3 M 30/i _1 Mpc runs at four epochs. Halos with 
Vmax > 200km/s were used to compute £,hh(r,z) for the 
SCDM and the rCDM models. A lower V max > 120km/s 
limit was used for the OCDM and ACDM 30 models. The 
difference in the V max limits is explained by the differ- 
As was shown, in the previous section, our £dm agrees very well with 
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Fig. 9. — The evolution of the correlation function of the dark 
matter (solid lines) and halos (dot-dashed line) for the SCDM 
model in the 30/i _1 Mpc AP 3 M simulation. The panels show 
the correlation functions at different redshifts. Only halos with 
maximum circular velocity larger than 200km/s were used to 
compute the halo correlation function. The number of halos 
used to estimate the correlation function is indicated in each 
panel. The poisson errors (dotted) and bootstrap errors (solid) 
are shown by vertical bars (see § 5.3 for details). 
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Fig. 10. — The same as Figure 9, but for the rCDMmodel. 
As in Figure 9, only halos with V ma x > 200km/s were used. 



ence in the matter density Oo that results in different 
mass resolution of the simulations. In the OCDM and 
the ACDM 30 runs halos with V max > 200km/s are scarce, 
while the poorer mass resolution does not allow us to re- 
duce the limit to 120km/s in the SCDM or the rCDM 
runs. The total number of halos found in the simulations 
is indicated in each panel. The number depends on the 
epoch and model and varies from the maximum of 1942 
(in the SCDM run at z = 1.9) to the minimum of 611 
(in the rCDM run at z = 3.6). The statistics of halos 
are poorer than in the ACDMgo simulation. Therefore, we 
plot the errorbars associated with each point of £hh ■ We 
estimate both poisson and bootstrap errorbars and plot 
the largest of the two. The bootstrap error bars have been 
estimated as follows. For each run and each epoch we have 
drawn 5 randomly selected samples of halos from the cor- 
responding halo population. The number of halos in each 
sample is one half of the total. We then compute the rms 
fluctuation between the samples and divide it by y/2 to 
get the la errorbar. 

Figures 9 — 12 show that in all models the evolution of 
the correlation function is qualitatively similar to that ob- 
served in the ACDM model. For example, the shape of the 



£hh is similar (power-law) in all models and is always dif- 
ferent from the corresponding £<j TO . This means that the 
scale- dependent bias is universal in cosmological models. 
Details of the evolution arc, however, model dependent. 
The most drastic differences are seen at the highest red- 
shifts. The figures, for example, clearly show that the bias 
at z = 3.4 has very different values in different models. 
While bias in the two low-f2 models is very similar, the 
distribution of halos at this redshift is only weakly biased 
in the SCDM model, as opposed to the strongly biased 
distribution in the rCDM model. The evolution of the 
halo correlation function at three different comoving scales 
(0.3, 1, and 3/i _1 Mpc) for all models is plotted in Figure 
13. The solid lines in both panels represent results for 
the ACDM model with the ART code in the 60/i _1 Mpc 
box. In the upper panel, the evolution is shown for the 
halo catalogs with a fixed number density of halos, which 
was achieved by varying the V m ax limit in different mod- 
els. Note that in this case we compare correlations of halo 
samples with different mass functions: the ACDM and 
OCDM halo samples contain many low-mass halos, while 
samples in the Oo = 1 models contain only massive ha- 
los. Such comparison is interesting for comparisons with 
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Fig. 11.— The same as Figure 9, but for the ACDM model. 
Only halos with Vmax > 120km/s were used to compute the 
correlation function (see § 5.3). 

observations when we know the number density of objects 
in the sample rather than their mass (or type). At scales 
<; l/i _1 Mpc, differences between the models are not sig- 
nificant. At 3/i _1 Mpc and at high z the amplitude of £hh 
in the SCDM model is significantly lower than in other 
models. The amplitude in the rest of the models is sur- 
prisingly similar. Therefore, if the biased galaxy formation 
scenario is correct and galaxies can be associated with host 
halos, this result may have interesting implications for the 
interpretation of clustering observations. To be able to dif- 
ferentiate between the models, we must know what type 
of the objects was used to estimate the clustering signal. 
The knowledge of the number density of objects in the 
sample is not sufficient. The point is to some extent illus- 
trated by the lower panel of Figure 13, where we compare 
the evolution of the £hh amplitude in ACDMgo, SCDM, 
and rCDM models for the halo catalogs with the same 
selection criterion (V max > 200km/s). It is obvious that 
in this case the differences between the models are signifi- 
cant. Although the differences are smaller at low redshifts, 
at z 3.5 the difference in the amplitude between ACDM 
and SCDM models is almost an order of magnitude. This 
difference can probably be explained by the delayed forma- 
tion of galaxy-size halos in the ACDM model as compared 
with the SCDM model. The halos in the ACDM form at 
lower redshifts with high statistical bias, while halos in the 
SCDM form systematically earlier and thus have had time 
to go through merging evolution. The effect of the latter 
is to decrease the bias (e.g., Moscardini et al. 1998). Note 
also that merging rates are higher in Q = 1 models (e.g., 
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Fig. 12.— The same as Figure 9, but for the OCDM model. 
As in Figure 11, only halos with Vmax > 120km/s were used 
(see § 5.3). 



Carlberg 1990 ). These results show that predictions of cos- 
mological models are very different for samples of objects 
selected with the same set of criteria for all models. 

The evolution of bias at scales 0.3 and l/i -1 Mpc is shown 
in Figure 14 for all models. Here again we compute the 
halo correlation function for the fixed number density of 
halos. Evolution of bias in all models is qualitatively simi- 
lar to that of the ACDM model discussed above: the bias is 
a very strong function of redshift. However, unlike the t;hh 
amplitude, the value of bias at these scales is very differ- 
ent among the models. This is not very surprising because 
when the number density of halos is fixed, different mod- 
els have very similar amplitudes of £hh (see Fig. 13), but 
very different amplitudes of £d m . The latter is explained 
by the differences in the cosmological parameters, normal- 
ization, and the shape of the power spectrum. A more 
interesting implication of the Figure 14 is that differences 
in bias get smaller at low redshifts, virtually disappearing 
at z = 0. The same effect can be observed in the evolution 
of the amplitude in Figure 13. As we will argue in the next 
section, the evolution of the halo correlations and bias at 
these scales is likely to be driven by the halo dynamics 
within nonlinear structures, in which case the differences 
between different cosmologies are largely erased. The evo- 
lution shown in Figures 13 and 14 provides, therefore, in- 
direct support for this point: the differences in clustering 
amplitude and bias between the models disappear at z <J 1, 
where most of the clustering signal comes from the halos 
located in nonlinear structures. 



EVOLUTION OF BIAS 



17 



1000 



100 



£10 



1 r 



1000 



100 



£10 r 



1 7 



"holo= 021 1,3 M P C " 3 




r-3.0 Mpc/h 



□ 



V > 200 km/sec 



r=0.3 Mpc/h 



r=3.0 Mpc/h 



■ ACDM box= 60 Mpc/h 



1+z 

Fig. 13. — The evolution of the halo correlation function at 
various scales for all models. The correlation functions shown 
in the upper panel were computed using a fixed number density 
of halos (implying different limits on the maximum circular ve- 
locity cut). The correlation functions shown in the lower panel 
were computed using all halos with V m ax larger than 200km/s. 
Results for only two 30 /i _1 Mpc simulations are presented in 
this panel. The other two simulations had too small numbers 
of halos. The solid lines in both panels represent results for the 
ACDMmodel with the ART code in the 60/i _1 Mpc box. 
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Fig. 14. — The evolution of bias for different models at two 
scales, r = l.0h~ Mpc (upper panel) and r = 0.3/i _1 Mpc 
(lower panel). The markers show results with the AP 3 M code, 
the solid curves represent results for the ACDM model with the 
ART code in the 60/i _1 Mpc box. 



6. DISCUSSION 

It is interesting to compare the evolution of the halo cor- 
relation function and the bias observed in our simulations 
with predictions of the analytical models and results of 



previ ous numerical siiuulaLiuiis. — The fact, LhaL ctttstermg 
stren gth of haloo at high rodohifto ia comparable to that 
at tha prooont e poch, hao been no ted in rcoulto of many 



simul ations (e . g . , Davis et al, 1985 1 Brainord fc Villumscn 



1994; Coh'n et al. 1997; and references therein). Bagla 
(1998) summarizes the generic behaviour of the correla- 
tion amplitude of halos above certain mass: the ampli- 
tude is high at very high redshifts, when halos are being 
formed, decreases thereafter and reaches a minimum, and 
then increases slowly and steadily until the present epoch. 
The results presented in the previous section (see Figs. 6 
and 13) are in agreement with this picture. Thus, there 
seems to be a good qualitative (although, in some cases, 
not quantitative) agreement among results of different nu- 
merical simulations concerning the evolution of the halo 



correlation function. 

Analytical models have reached a sufficient degree of so- 
phistication to be able to predict the evolution of halo clus- 
tering in mildly nonlinear regimes (see § 2). The hal os are 
found to for m at the peaks of the density field (e.g., Frenk 



et al. 1988 ) and their bias exhibits a simple sca ling rela- 
tion with t he height of these peaks (Kaiser 1984; Bardeen 
et al. 1986| ). At any given epoch the halo population rep- 
resents a mix of halos formed at different redshifts: newly- 
born or already evolved through merging. The evolution 
of the correlation of halos in such a hierarchical framework 
is described using the extended Press-Schechter formalism 
(MW). To compare predictions of the analytical models 
with our results, we will use the approximation to the evo- 
lution of the effective bias given by eq. (4) (Moscardini et 
al. 1998). This approximation describes the evolution of 
bias of a sample of all halos above a certain mass. This is 
roughly equivalent to our definition of halo samples with 
a limiting maximum circular velocity. The prediction of 
this approximation is shown in Figure 8 with the dotted 
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line, where we used 6 e //(0) = 0.51 and (3 — 1.90 (see 
eq.[4]) appropriate for our ACDM model and for the mass 
limit of M > 10 11 /i~ 1 M Q (Moscardini et al. 1998). The 
analytical model is expected to provide a good approxi- 
mation at scales where £hh & 1 (Mo et al. 1996), i.e. at 
r ;> 4 — 5/i _1 kpc (see Table 4). Therefore, the analyti- 
cal prediction should be compared with the curve showing 
evolution of bias at r = 4.8/i -1 Mpc. 

The comparison shows that both the numerical result 
and the analytical model predict a rapid decrease of bias 
with cosmic time. Moreover, at high redshifts (z ;> 3) they 
agree well quantitatively. At lower redshifts, however, the 
two predictions deviate from each other and are different 
by a factor of 2 at z = 0: at r = 4.8/i _1 Mpc almost no 
bias is observed in the simulation, while strong anti-bias 
of b « 0.5 is predicted by the analytical model. The differ- 
ences are not surprising, given the differences in our defi- 
nition of a halo from that of the Press-Schcchter halo. The 
definition of the latter does not include "satellite" halos; a 
halo ceases to exist, once it becomes bound to another halo 
(i.e., "merges") and orbits inside that halo's virial radius. 
Our definition, on the other hand, does take satellite ha- 
los into consideration, because we include in our halo list 
every gravitationally bound clump of particles, regardless 
of whether it is also bound to a larger system or not. The 
similarity between our numerical result and the prediction 
of the analytical model at z J> 3 is then an indication 
that the two definitions are equivalent at these high red- 
shifts. Indeed, large systems such as clusters and groups 
have not yet formed at these redshifts, and the fraction of 
satellite halos in our catalogs (i.e., all halos above a cer- 
tain mass limit) is relatively small. At smaller redshifts 
the ever larger fraction of halos become satellites to more 
massive halos and the two halo definitions result in rather 
different halo samples. This explains the large difference 
predicted for the value of bias at z = 0. We have found 
that in our simulation the z = amplitude of the 2-point 
correlation function of the Press-Schecter halos (i.e., iso- 
lated in terms of their virial radius) at r = 4.8/i _1 Mpc is 
approximately two times as small as the amplitude of the 
correlation function of the BDM halos shown in Fig. 7, 
resulting thus in the twice as small bias of b « 0.5. In this 
respect, we believe that there is no contradiction between 
our results and predictions of the analytical model, once 
the difference in the halo definition is taken into account. 
This also explains the difference in the bias value with the 
result of Jing (1998), who found bias of b « 0.5 — 0.7 for 
the galaxy-size halos defined with the FOF algorithm 17 . 
We find also that, contrary to the assumption of the an- 
alytical models, the small-scale bias of the halo distribu- 
tion is scale- dependent regardless of the halo definition. 
The large-scale bias (r ^> 5/i _1 Mpc), on the other hand, is 
not probed in our simulations and may be independent of 
scale, in accord with assumptions of the analytical mod- 
els. This indeed was suggested by the simulations of Jing 
(1998). We think that the most encouraging result is the 
agreement between numerical and analytical modeling on 
the general form of the bias evolution demonstrated in 



Fig. 8. This indicates that we now have a solid general 
understanding of the nature of bias and of the processes 
driving its evolution at redshifts z J> 2 — 3. 

At smaller redshifts, the merging rate is considerably 
smaller and small-scale correlations are sensitive to the 
dynamics of halos inside the nonlinear structures. Partic- 
ularly, the dynamics and the clustering evolution of satel- 
lite halos in high-density regions are essentially indepen- 
dent of the background cosmology and are driven by such 
processes as dynamical friction and tidal stripping (e.g., 
KGKK). These processes tend to suppress the growth of 
the correlation amplitude, thus counteracting the cluster- 
ing growth due to the gravitational pull. This leads to the 
anti-bias observed in our simulations at nonlinear scales 
and at small redshifts. Indeed, the correlation amplitude 
at small scales (r <J 5ft. _1 Mpc) is approximately constant 
at redshifts z ss — 1 (see Figs. 6 and 13). Moreover, 
anti-bias is observed in all cosmological models studied in 
this paper (see Figs. 9-12). The fact that differences in the 
correlation amplitude and bias, existing between the cos- 
mological models at z ~ 3, virtually vanish by the present 
epoch (see Figs. 13 and 14) argues for the importance of 
the nonlinear halo dynamics. This result also implies that 
low-redshift clustering depends only weakly on the back- 
ground cosmology. Therefore, the information about the 
underlying cosmological model can probably be extracted 
only from the high-z (z J> 3) clustering data. As was dis- 
cussed in § 4.3, the numerical resolution required to assure 
the survival of halos in high-density regions is high, and 
has not been reached in previous simulations. Our results 
therefore indicate that high resolution is important for cor- 
rect modelling of the bias evolution at small redshifts. 

Although our definition of a halo is different from that 
used in the conventional Press-Schechter framework, we 
believe that it is closer to what can be identified in a 
simulation as a galaxy location. It seems likely that in 
every sufficiently massive (M £ lO n /i _1 M ) gravitation- 
ally bound halo baryons will cool, form stars, and produce 
an object ressembling a galaxy (e.g., Kauffman, Nusser 
& St cinmctz 1997; Roukcm a et al. 1997; Yepes et al. 
1997; ISalucci k Persic 1997| ). We believe, therefore, that 



each of the halos in our catalog can be associated with a 
"galaxy" 18 . Observationally, many distinct galaxies are 
located well inside the virial radii of massive galaxies, 
groups, and clusters. These galaxies, non-existent by def- 
inition in the "virial overdensity" halo catalogs, are in- 
cluded in galaxy surveys and are used to compute the 
correlation function. Our definition, therefore, is natu- 
ral, if the goal is to compare the observations with predic- 
tions of the numerical simulations. Although, as we dis- 
cussed above (§ 4.3), the resolution (and, correspondingly, 
the computational costs) required for galaxy-size halos to 
survive in the tidal fields of high-density regions is quite 
high, the subsequent comparison with the observations is 
straightforward and does not require ambiguous correc- 
tions for the "overmerging" . 

In §5.2 we have compared the correlation function of 
halos S,hh in our 60/i _1 Mpc ACDM simulation with the 



17 In this case again the halos are identified as isolated density peaks with an overdensity exceeding the value expected of the virialized object. 
No satellite halos can be identified with this definition. 

18 We cannot, of course, unambiguously assign type, color, luminosity, or other galactic properties to our halos without additional modelling. 
Our results, therefore, should be applicable to "global" galaxy surveys, such as the APM survey, in which galaxies are selected solely on the 
basis of their luminosity (loosely related to the mass). 
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z = correlation function of galaxies £ gg (Baugh 1996) in 
the APM catalog (Loveday et al. 1995) and found a very 
good agreement between the two. The APM galaxy corre- 
lation function is measured very accurately, which makes 
the agreement within the lcr-errorbars very striking (see 
Fig. 7). Recently, the advent of new faint galaxy surveys 
allowed the measurement of the clustering surveys at red- 
shifts of z m — 1 (e.g., Le Fevre et al. 1996; Shepherd et 
al. 1997; Carlberg et al. 1997; Connolly et al. 1998; Post- 
man et al. 1998; Carlberg et al. 1998). Unfortunately, 
there is some disagreement between these studies concern- 
ing the amplitude of the clustering signal, which possi- 
bly indicates that there is morphology and/or luminosity 
segregation in the clustering of the intermediate redshift 
galaxies. It is indeed true that the correlation amplitude 
depends on the luminosity of the galaxies (e.g., Carlberg 
et al. 1998; Postman et al. 1998). The comoving correla- 
tion length ro measured in these surveys varies, generally 
lying in the range tq <~ 2 — 4/i~ 1 Mpc with an average value 
of about ro ~ 3/i _1 Mpc. The correlation length of bright 
galaxies is, however, somewhat larger and consistent with 
the correlation length of the local galaxies ro ~ 5/i _1 Mpc 
(Carlberg et al. 1998; Postman et al. 1998). Comparison 
of the above values of ro with correlations in our ACDM 
simulation (Table 4), shows that we predict a correlation 
length in good agreement with that of bright galaxies and 
somewhat larger than the r value of the faint galaxies. 
The latter, however, is measured with rather large uncer- 
tainties and our values of ro are actually consistent with 
the observed ones within lcr. 

Remarkable progress of the high-redshift galaxy detec- 
tion techniques, based on the search for the signatures 
of the Lyman break in the colors of faint galaxies (Stei- 
del & Hamilton 1992, 1993), resulted in a rapid growth of 
amount and quality of the clustering observations at z m 3 
(Steidel et al. 1998; Giavalisco et al. 1998; Adelberger et 
al. 1998). These "Lyman break galaxies" (LBGs) were 
found to be clustered at z ~ 3 as strongly as the present 
day galaxies. The real-space correlation function of these 
galaxies was well-described by the conventional power-law 
form with the value of the slope 7 and the correlation 
length r consistent with the z s=a values (Giavalisco et 
al. 1998): 7 w 1.7 - 2.2 and r w 3/i _1 Mpc. These values 
are in reasonably good agreement with z = 3 values for 
our V max > 120km/s catalog (Table 4). The value of bias 
that is measured for the LBGs at scales of ss 5 — 10/i _1 Mpc 
is b ~ 1.5,3.6,4.5, for the fl = 1, Clo = 0.2 (open), and 
Oo = 0.3 (flat) models, respectively (Giavalisco ct al. 1998; 
Adelberger et al. 1998). These values can be compared 
with b(z = 3) for halos in our simulations, shown in Figs. 8 
and 14. All models agree with the observations, within the 
uncertainties of the galaxy-to-halo mapping. This result 
is in general agreement with the results of the numerous 
recent numerical studies that modeled the clustering of 
LBGs (e.g., Wechsler et al. 1998; Jing & Suto 1998; Bagla 
1998; Governato et al. 1998). However, there appear to 
be some puzzling details in comparisons with the data. 
The observed value of &lbg ~ 3 — 4 can be reproduced 
in the ACDM for massive halos with V max > 200km/s 
(see top panel of Fig. 8). This is in agreement with al- 
most all other theoretical studies. However, the corre- 
lation function of the LBGs measured by Giavalisco et 
al. (1998), does not agree with the correlation function 



measured for the V max > 200km/s of our ACDM ha- 
los: ro ~ 3ft. _1 Mpc for LBGs versus ro s=s 4.5/i _1 Mpc for 
the halos. This disagreement was actually noticed in the 
study of Adelberger et al. (1998) who used the count- 
in-cells analysis to derive the value of bias. The values 
of the parameters of the correlation function that were 
derived from the observed rms fluctuations of galaxies in 
cells of w 12/i _1 Mpc {(Tgai ~ 1.1 ± 0.2), are considerably 
higher than those measured directly by Giavalisco et al. 
(1998). The corresponding value of ahaioir = 12/i _1 Mpc) 
in our ACDM simulation is 0.6,0.7,0.9 for halos with 
V max > 120, 150, 200km/s, respectively. This is consistent 
with the interpretation of the LBGs as objects residing in- 
side massive (V max > 200km/s or M £ lO 12 /?.-^©) DM 
halos. This result supports the interpretation of Adel- 
berger et al. (1998) and suggests that the correlation am- 
plitude of the LBGs may be higher than that obtained 
from the observed angular correlation function (Giavalisco 
et al. 1998). 

Overall, we believe that the comparisons discussed 
above indicate that there is good agreement between our 
results and the clustering data at both low and high red- 
shifts. This implies that hierarchical models in which ob- 
served galaxies form in the host DM halos naturally ex- 
plain the observed galaxy clustering at different epochs, 
including excellent agreement with the accurately mea- 
sured z = correlation function. On the other hand, 
the generic form of the bias evolution observed in the nu- 
merical simulations at high redshifts agrees well with the 
prediction of the analytical models based on the extended 
Press-Schechter formalism. This implies that we under- 
stand the nature of the bias and the processes that drive 
its evolution at high z. At low redshifts, the bias evolution 
of gravitationally bound halos is driven by the dynamical 
processes inside the nonlinear structures which are largely 
independent of cosmology. The study of these processes is 
important for a successful modelling of galaxy clustering 
at z & 1. 

7. SUMMARY 

We have studied the evolution of the correlation func- 
tion and bias of galaxy-size halos in different cosmological 
models (ACDM, OCDM, rCDM, and SCDM). The high- 
resolution of our numerical simulations allowed us to avoid 
the overmerging in the high-density regions and estimate 
the correlation amplitude and bias directly at small (down 
to ~ 100/i _1 kpc) scales. The main results and conclusions 
presented in this paper are as follows. 

1. At all epochs, the 2-point correlation function of 
galaxy-size halos ^hh is well approximated by a power-law 
£hh = (f/fo) 7 w ith the slope 7 w 1.6 — 1.8. The corre- 
lation length r at z — is w 5/i _1 Mpc, regardless of the 
minimal mass limit of the halo samples. At high redshifts, 
the correlation function evolves non-monotonically: ro de- 
creases somewhat between redshifts of 5 and 3, and then 
increases steadily until z = 0. For the most massive halos, 
the correlation length at z w 5 is comparable to that at 

2 = 0. 

2. The difference between the shape of the £hh and the 
shape of the correlation function of matter results in a 
scale-dependent bias at scales <; 7/i~ 1 Mpc. Wc find this 
to be a generic prediction of the hierarchical models inde- 
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pendent of the epoch and of the model details. 

3. Another generic prediction is that the comoving am- 
plitude of the correlation function for halos above a certain 
mass evolves non-monotonically: it decreases from an ini- 
tially high value at z ~ 3 — 7, and very slowly increases at 
z <; 1. This behaviour at large scales was demonstrated 
by a number of authors (see § 6). Here, we have shown 
that this behaviour also applies to the correlation ampli- 
tude at small scales (< lft," 1 Mpc). The non-monotonic 
evolution of the correlation function calls into question 
the usefulness of the simplistic "e-models" as a descrip- 
tion of the clustering evolution. We note, however, that at 
z <j 1 the evolution of the halo correlation function is ap- 
proximately monotonic (albeit dependent on scale). The 
very slow evolution of the halo correlation amplitude in 
comoving coordinates at these redshifts implies a value of 
e i=a — 1, which is in agreement with the values preferred 
by the observations (e.g., Postman et al. 1998). 

4. The evolution of the halo correlation function is very 
mild compared to evolution of the dark matter correlation 
function. The latter evolves by a factor of ~ 10 — 60 (de- 
pending on scale) between redshifts of w 7 and 0, while 
the difference in amplitude of the former between any two 
epochs is less than a factor of 2. The large difference in 
the evolution rates of the matter and halo correlation func- 
tions means that the bias evolves rapidly with cosmic time: 
it changes from high b values of ^ 2 — 5 at z ^ 3 — 7 to 
anti-bias b values of ~ 0.5 — 1 on small <; 5/i _1 Mpc scales 
at z = 0. 

5. We find that our results agree well with existing 
clustering data at different redshifts, indicating a general 
success of the hierarchical models of structure formation 
in which galaxies form inside the host DM halos. Particu- 
larly, we find excellent agreement in both slope and ampli- 
tude between £,hh{z = 0) in our ACDMgo simulation and 
the galaxy correlation function measured using the APM 
galaxy survey. At high redshifts, all models reproduce well 
the observed clustering of the Lyman-break galaxies. Our 
results imply that for high-redshift clustering to be used 
as a cosmological test, it is crucial that we know what type 
of objects are used to estimate the clustering signal. The 
knowledge of the number density of objects is not suffi- 
cient (see § 5). 

6. We find good agreement at z p> 2 between our results 
and predictions of the analytical models of bias evolution 
(MW; Matarrese et al. 1997). This indicates that we now 
have a solid understanding of the nature of the bias and 



of the processes that drive its evolution at these z. We ar- 
gue, however, that at lower redshifts the evolution of bias 
is driven by dynamical processes, i.e., dynamical friction 
and tidal stripping, inside the nonlinear high-density re- 
gions such as galaxy clusters and groups. These processes 
do not depend on cosmology and tend to erase the differ- 
ences in clustering properties of halos that exist between 
cosmological models at high z. The latter result implies 
that low-redshift clustering is probably not a very strong 
discriminator between cosmological models. 

We believe that the success of the current theoretical 
models in interpreting the clustering data forms a solid 
foundation for further sophistication of the models by in- 
cluding the processes important for galaxy formation (such 
as dynamics of baryons, cooling, star formation, and stel- 
lar feedback). These models would allow one to predict 
the observed properties of galaxies and thus mimick the 
observational selection criteria, allowing for a robust com- 
parison between the model and the data. We believe 
that the differences between high-z clustering properties 
of objects in different cosmological models demonstrated 
in this study (see § 5), the improved theoretical models, 
and the ever increasing amount of new clustering data can 
be successfully combined in the near future to put useful 
constraints on the cosmological parameters describing our 
Universe. 
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